Mathematical Structure of Evolutionary Theory 
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Here we postulate three laws which form a mathematical framework for the evolutionary dynamics 
in biology. The second law is most quantitative and is explicitly expressed in a unique form of 
stochastic differential equation. Salient features of Darwinian evolutionary dynamics are captured 
by this law: the probabilistic nature of evolution, ascendancy, and the adaptive landscape. Four 
dynamical elements are involved in the present formulation: the ascendant matrix, the transverse 
matrix, the Wright fitness function, and the stochastic drive. The first law may be regarded as a 
special case of the second law. It gives the reference point to discuss the evolutionary dynamics. The 
third law describes the relationship between the focused level of description to its lower and higher 
ones, and also defines the dichotomy of deterministic and stochastic drives. A precise description of 
Wright's adaptive landscape is given and a new interpretation of Fisher's fundamental theorem of 
natural selection is provided. The generality of the proposed laws is supported by the demonstration 
of their equivalence to a general non-equilibrium dynamics formulation, based on a recently proved 
theorem. Though the proposed laws are not the rigorous description, they should be regarded 
as mathematically what the laws of evolution might be. In addition, they provide a coherence 
framework to discuss several current evolutionary problems, such as speciation and stability, readily 
showing that statistical physics tools can be applied to Darwinian dynamics. 



I. INTRODUCTION 

Progresses in experimental biology and new data from 
field observation after the neo-Darwinian synthesis pose 
new questions to be answered and call the attention to 
previously unanswered questions. Biologists have been 
responding to this demand with tremendous activities, 
ranging from a critical discussing 1 of existing theories 
and courageous exploring 2 additional one, to the serious 
consideration of the epistasis of gene interactions 3 , to a 
reevaluation of shifting balance process 4-6 , to a reexam- 
ining the concept of species 7 and of the fundamental theo- 
rem of natural selection 8,9 , to various elegant mathemat- 
ical models of speciation 10-12 . Even new philosophical 
implications were speculated 13 . Such efforts have always 
enriched the theoretical and conceptual understanding of 
evolution. 

Among the various approaches, it has been noted 3 that 
a quantitative formulation of evolution dynamics may be 
of more importance to answer the new questions. Ver- 
bal description has been found to be inadequate, because 
many contributions to evolutionary dynamics, either in- 
dependent or interactive, are of same magnitude, and 
the interactions can be highly nonlinear. It was from 
this consideration Stewart 10 built his model based on 
symmetry-breaking and Gavrilets 11 advanced his holey 
adaptive landscape model. This line of reasoning is fur- 
ther developed in the present article. Here we make an 
attempt to formulate a general and quantitative mathe- 
matical framework which appears broad enough to incor- 
porate the ideas of Stewart 10 and Gavrilets 11 and is con- 
ceptually consistent with the Darwinian dynamics. Two 
of the most influential concepts in evolutionary biology, 
the adaptive landscape and the fundamental theoretical 
of natural selection, are built naturally into the present 



formulation. The present work may also be regarded as 
an attempt to unify approaches from both biological and 
physical sciences. 

The present mathematical approach is based on the 
continuous approximation which treats the populations 
as continuous variable. This approximation has been well 
studied in biology, documented, for example, in both a 
commentated collection of historical articles 14 , in a recent 
monograph 15 and in a recent textbook 16 , and in an online 
book 17 , and has been successfully employed in population 
genetics. In present article we will not elaborate further 
on this approximation. This implies that the equations to 
be discussed are of differential equation type. To be more 
precise, we will postulate three laws for evolution and the 
most important law, the second law, will be expressed in 
a unique form of stochastic differential equation. The 
connection of the present triad laws to previous results 
will be discussed. 

We organize the rest of article as follows. In section 
II we postulate and discuss the three laws of evolution. 
Four dynamical elements are needed in this formulation. 
In section III the connection of the postulated three laws 
to previous formulations is discussed. An explicit compu- 
tation of this connection is demonstrated in section IV. 
We also discuss the compatibility of present formulation 
with two current speciation models. The implications of 
present formulation were discussed in section V, and we 
conclude in section VI. 



II. LAWS OF EVOLUTION 

In this section we postulate and discuss three laws on 
evolutionary dynamics. They form a quantitative math- 
ematical framework for evolutionary dynamics. Then we 
discuss Fisher's fundamental theorem of natural selection 
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and conclude that it is indeed an indispensable relation- 
ship. We start with the most important and quantitative 
law, the second law. 

A. Second Law 

The central question naturally arises that how do we 
describe the evolutionary dynamics quantitatively and 
what are the dynamical elements? To be specific, let 
us consider an n component biological system. The n 
components may be the species in a evolutionary game 18 , 
or the traits to describe the speciation 10 , or genes in the 
description 12 , or any quantities required to specify the 
system. The value of j component is denoted by qj. 
The n dimensional vector q T = (qi, Q2, q n ) is the state 
variable of the system. Here the superscript r denotes the 
transpose. The dynamics of state variable is described by 
its speed q t = d(\ t /dt moving in the state space. 

We postulate that the dynamics of the system is gov- 
erned by following special form of stochastic differential 
equation, which consists of four dynamical elements, the 
semi-positive definite symmetric ascendant matrix A, the 
anti-symmetric transverse matrix T, the scalar function 
called Wright fitness function ip, and the stochastic drive 

[A{d U t)+T{d t ,t)]q_ t = V^(q,,t)+^(q,,t) , (1) 

and supplemented by the following relationship: 

(£(q*,t)r(q*',*')> =2A(d t ,t)e5(t-t') , (2) 

and (£(qt, t)) = 0. The connection of these two equations 
to conventional approaches will be discussed in next sec- 
tion. To ensure the independence of the dynamics of 
each component, we assume det[A(q, t) + T(q, t)} ^ 0. 
Here the subscript t denotes that the state variable is a 
function of time, and 8{t) is the Dirac delta function. In 
Eq.(2) we have assumed the stochastic drive is Gaussian 
white noise with zero mean. Factor 2 is a convention 
choice for the present formulation, and e is a positive nu- 
merical constant, which for many situations might be set 
to be unity, e = 1, without affecting the biological de- 
scription. The relationship between the stochastic drive 
and the ascendant matrix expressed by Eq.(2) guarantees 
that the ascendant A(q t ,t) is semi-positive definite and 
symmetric. The average (...) is carried over the dynamics 
of the stochastic drive, and V is the gradient operator in 
the state space. 

It is straightforward to verify that the symmetric ma- 
trix term is 'ascendant': q[A(q t ,t)q t > 0. Its dynamical 
effect is to increase the fitness in terms of the Wright fit- 
ness function ip(qt, t). Here we point out that the Wright 
fitness function defines the adaptive landscape. Hence as- 
cendant matrix enables the systems to have the tendency 
to seek the largest fitness. This feature will be explicitly 
manifested in the discussion after the first law. The anti- 
symmetric matrix permits 'no change': qJ"T(qt,i)q t = 0, 



therefore conservative. Dynamically it will not change 
the fitness. A manifestation of the transverse dynam- 
ics is the oscillatory behavior. The dynamical effect of 
the stochastic drive £(qt,i) on fitness is random: It may 
either increase or decrease the fitness. With above inter- 
pretation, the static effect natural selection is represented 
by the gradient of fitness function, V^p{(\ tl t). The clear 
and graphical discussion of such Wright fitness function 
was first clearly expressed by Wright 19 in his concept of 
adaptive landscape, with which the present fitness func- 
tion is identified. The tempo of natural selection is repre- 
sented by the ascendant and transverse matrices. Eq.(l) 
states that the four dynamical elements, the gradient of 
fitness, the stochastic drive, the ascendant dynamics, and 
the transverse dynamics, must be balanced to generate 
the evolution dynamics. 

Eq.(l) is the fundamental equation of evolutionary 
dynamics expressed in a unique form of stochastic dif- 
ferential equation. In accordance with above discus- 
sion on stochastic drive and ascendant matrix, we may 
call the supplementary equation, Eq.(2), stochasticity- 
ascendancy relation, and will discuss it in the last sub- 
section in this section in relation to Fisher's fundamental 
theorem of natural selection. 

The Wright fitness function ip is similar to a potential: 
It is in fact opposite in sign to the typical potential en- 
ergy used in physical sciences. It is in general a nonlinear 
function of state variable. If it is further independent of 
time and is bounded above, the stationary distribution 
function p(q, t — oo) for the state variable, the proba- 
bility density to find the system at q in state space, is 
expected to be a Boltzmann-Gibbs distribution: 

/9 (q,i = oo) = |exp|^| , (3) 

with Z — J n™=i dqi exp {^(q)/e} the partition function, 
the integration over whole state space serving as the nor- 
malization factor. Its justification will be given in next 
section. It is interesting to note that the dynamical as- 
pects of evolution, the transverse and the ascendant ma- 
trices, do not explicitly show up in Eq.(3). From this ex- 
pression we observe that the larger the constant e is, the 
wider the equilibrium distribution would be, and more 
variation would be, or, the smaller the e is, the narrower 
the distribution. In this sense we may call e the evolution 
hotness constant. The existence of such a Boltzmann- 
Gibbs type distribution suggests a global optimization. 

There are a few immediate and interesting conclusions 
to be drawn here. Near a fitness peak, say at q = q pea fe, 
we may expand the Wright fitness function, V(q) = 
V'(qpeafe)- (q-qp e afe) T ^(q-qp e afc)/2 + 0(|q-q pe afe| 3 )- 
Here U is a positive definite symmetric matrix as a con- 
sequence at the fitness peak. The stationary probability 
density to find the system near this peak is of a typical 
Gaussian distribution: 
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Thus, away from the fitness peak, the probability to find 
the system will be exponentially small. Such a behavior 
has long been observed in many biological models 15 ' 17 . 

One may then wonder about how does the system move 
from one fitness peak to another? This process was first 
visualized by Wright 19 . The relevant mathematical cal- 
culation seems to be first done by Kramers 20 . It had 
been applied to biology 21 , where it was shown that the 
stochastic drive must be involved. Quantitatively, the 
hopping from one peak to another must be aided by the 
stochastic drive. The dominant factor in the hopping 
rate T is the difference in fitness between the peak and 
the highest point (saddle point q sa ddie) to cross the valley 
to another peak 20-22 : 



r cx exp < — 



V'(qpeafe) - IpiQsaddle) 



(5) 



This rate can easily be exponentially small. It is a quanti- 
tative measure of robustness and stability. Hence it may 
explain the usual observation, for example, that species 
is rather stable if viewing the peak as a definition for 
species. Nevertheless, Eq.(5) grants the possibility to 
hop between peaks when the stochastic drive is finite. 

The second law as expressed by Eq.(l) and (2) appar- 
ently capture the major features of the evolution dynam- 
ics first described by Darwin and Wallace 23 , exhaustively 
exposed by Darwin 24 , and further developed by Fisher, 
Haldanc, Muller, and Wright 25 , and by many others 15 ' 17 . 
Obviously it expresses the evolutionary process as a tin- 
kering process 26 . The necessity and chance 27 are repre- 
sented by the Wright fitness function and the stochastic 
drive, respectively. 



B. First Law 

The first law is a statement for the situation that there 
is no stochastic drive in the evolution, i.e., e = 0, the 
smallest possible value of the evolution hotness constant. 
This is clearly an approximation, because the variation 
is always there 28 , though it may be regarded to be small 
and slow under certain time scale. Allowing the stochas- 
tic drive be negligible, Eq.(l) becomes 



[A{q t ,t)+T{q t ,t)]i k = Vil>{q t ,t) . 



(6) 



Because of the ascendant matrix A is non-negative, the 
system will approach the nearby attractor determined by 
its initial condition, and stay there forever. Specifically, 
because qfiAfat, i)q t > and q[T(q t , t)q t = 0, Eq.(6) 
leads to 



q t - VV>(q*,t) > . 



(7) 



This equation implies that the deterministic dynamics 
cannot decrease the fitness: The speed of state variable 
q t is in the same direction of the gradient of the Wright 
fitness function V^(qt, t). If the ascendant matrix is pos- 
itive definite, i.e. qj"^4(q t , t)q t > for any nonzero q t , 
the fitness of the system always increases. Hence, the 
first law clearly states that the system has the ability to 
find the local fitness or adaptive peak determined by the 
initial condition. 

We have two remarks here. First, from the mathemat- 
ical theory of dynamical systems, there are in general 
three types of attractors 29 : point, periodic, and chaotic 
(strange). The point attractors have been well explored 
in evolutionary biology since the work of Wright, corre- 
sponding the fitness peaks. Other two types of attrac- 
tors have also been observed in biology 30 ' 31 . Second, if 
we further assume the ascendant matrix A = 0, the evo- 
lutionary dynamics will not change the system's fitness, 
hence is completely conserved. This is precisely what 
can be obtained from the Newtonian dynamics 32 . Based 
on this consideration one may incline to conclude that 
Newtonian dynamics is a special case of the Darwinian 
dynamics as expressed by Eq.(l) and (2), a conclusion 
many biologists may consider obvious 33 ' 34 . 

The tendency implied in Eq.(6) to approach an at- 
tractor and to remain there has been amply discussed by 
Aristotle, well known in physics and biology. It is evident 
that in the present setting the first law is mathemati- 
cally a special case of the second law. Nevertheless, this 
law does give us a reference point to define species and 
other relevant quantities in a clean manner, if stochastic- 
ity could be ignored. 



C. Third Law 

The third law is a relationship law. It allows us to 
define the connection of the current level of description 
to its lower and higher ones. It is a reflection of the hi- 
erarchical structure of the whole dynamics. The essence 
of this law is to acknowledge the existence of two time 
scales: the macro time scale with which we observe the 
system dynamics at the focused level of description and 
the micro time scale with which fine structures and lower 
level dynamics come into play. 

Specifically, it may be stated as follows: The Wright 
fitness function ip(q_, t) has the contribution from lower 
level in terms of time average on the time scale of cur- 
rent level, the contribution from the interaction among 
various components of the current level, and the contri- 
bution from higher level. The stochastic drive £ (q, t) is 
the remainder of all those contributions whose dynamics 
is fast on the time scale of current focus. Hence its av- 
erage in time is zero. This stochastic contribution may 
be either unknown from a more fundamental level or un- 
necessary to be specified in details. Only its probability 
distribution is needed and is approximated by a Gaus- 
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sian distribution in the present article. The stochastic 
drive determines the ascendant matrix A(q, t), and the 
transverse matrix T(q, t) should be further determined 
by the dynamics of the system. 

The lower level contribution to the Wright fitness func- 
tion ip(q, t) and the stochastic drive £(q, t) may allow us 
to compute the intrinsic fitness landscape and the intrin- 
sic source of evolution. However, historically the compu- 
tation of this contribution tends to neglect the horizontal 
interaction among different components, which is usu- 
ally nonlinear. On the other hand, the same and higher 
level contributions may suggest that a control mecha- 
nism, such as a feedback, may be from both of them 
in a large perspective. The combination of all three of 
them suggests that the evolution is nonlinear, asymmet- 
ric, mutually interactive, and stochastic, and may be con- 
trollable. 

There is a degree of uncertainty and arbitrariness in 
the assignment of different levels of descriptions and 
the dichotomy of deterministic and stochastic terms in 
Eq.(l). This dilemma has been amply discussed in both 
physical 22 and biological 35 sciences. This has also been 
reflected in the mathematical theory of stochastic pro- 
cess. Our way to solve this problem will be proposed in 
next section in connection to usual dynamics, which is 
different from both Ito and Stratonovich approaches and 
may be interpreted as an indication of the richness of the 
hierarchical structure. 



D. Fundamental Theorem of Evolution 

In his classic treatise on genetical foundation of 
evolution 36 , Fisher stated his fundamental theory of nat- 
ural selection: The rate of increase in fitness of any or- 
ganism at any time is equal to its genetic variance in 
fitness at that time. This is truly a remarkable insight: 
It immediately connects the adaptation in fitness land- 
scape to the stochastic drive. 

As implied in the first law, the ascendancy of the sys- 
tem is described by the ascendant matrix A, which in 
turn is completely determined by the stochastic drive ac- 
cording to the stochasticity-ascendancy relation, Eq.(2). 
The discussion followed Eq.(5) indicates that the ability 
of system to find a better fitness peak, not only the local 
fitness peak, or, to reach the global equilibrium, is guar- 
anteed by the stochastic drive. This suggests that Eq.(2) 
is a statement on the unification of the two completely 
opposite tendencies: adaptation and randomization. 

Ever since Fisher's proposal of the fundamental the- 
ory of natural selection, misrepresentation and misun- 
derstanding have been associated with this insightful 
statement, as discussed by Crow 8 and Grafcn 9 . Though 
Fisher might have confined his discussion within ge- 
netic variations, the comparison of Fisher's statement to 
Eq.(2) points to a remarkable similarity: On the left hand 
side of Eq.(2) is the measure of the variation, and on the 



right hand side may be interpreted as the rate to fitness 
peak with a proper choice of unit. Hence, we may also 
call Eq.(2) the fundamental theorem of evolution in a 
tribute to Fisher's great insight, or simply the Fisher's 
theorem. Its importance is evident: It is an integral part 
of the second law of evolution, the Darwinian dynamics. 
It enables the system to find the global maximum fitness 
peak. Nevertheless, it should be pointed out that the fun- 
damental theorem of evolution in the present formulation 
is independent of the Wright fitness function. 

Note: In Fisher's original formulation^ 6 he listed sev- 
eral exceptions to his fundamental theorem. According 
to the sharp analysis of S.J. Gould in his The Struc- 
ture of Evolutionary Theory (2002) Fisher's exceptions 
completely negate the importance of Fisher's theorem! 
The present author believes that Fisher's misinterpre- 
tation of his theorem is the source of long confusion 
and debate in population genetics continuing to these 
days. See for example a most recent text book 16 . In 
the present manuscript the precise mathematical formu- 
lation and clear ascendant matrix interpretation make all 
Fisher's exceptions unnecessary. Hence, the present for- 
mulation retains the fundamental nature of this theorem. 

The debate on Fisher's theorem is interesting even 
from physics point of view: There has been an enormous 
amount of confusion and debate on something similar in 
non- equilibrium thermodynamics. Perhaps the above for- 
mulation and the demonstration of the equivalence in the 
next section have been done in physics literature. Sugges- 
tions on relevant literature from readers would be highly 
appreciated. 

III. CONVENTIONAL FORMULATION 

A. Standard Stochastic Differential Equation 

Now we make the connection between the dynamics 
described by Eq.(l) and (2) to the dynamical equations 
typically encountered in evolution. We start with the 
standard stochastic differential equation: 

q* = f(q*,*) + C(q*,f)- (8) 

Here f (q, t) is the deterministic nonlinear drive of the sys- 
tem, which includes effects from both other components 
and itself, and the stochastic drive is £(q, t), which differs 
from that in Eq.(l) but is governed by the same dynam- 
ics. For simplicity we will assume that f is a smooth 
function whenever needed. Without loss of generality 
the stochastic drive in Eq.(8) is assumed to be Gaussian 
and white with the variance, 

(t(dt,t)C(ci t ,,t'))=2D(c lt ,t)e6(t-t'), (9) 

and zero mean, (£(q t ,t)) = 0. Eq.(9) is consistent with 
Eq.(2). Again here (...) indicates the average with respect 
to the dynamics of the stochastic drive. According to the 
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biology convention the semi-positive definite symmetric 
matrix D = {-D^} with i,j — 1,2, ...,n is the diffusion 
matrix. Both the divergence and the skew matrix of the 
nonlinear drive f are in general non-zero: 

V-f^O, Vxf^O. (10) 

Here the divergence is explicitly V • f = Y^j=i dfj/^Qj = 
tr(F), and the skew matrix V x f is twice the anti- 
symmetric part of the selection matrix S: (V x f)y = 
S 3 i - Sij with Sij = dfi/dqj , i,j = 1,2, ...,n. The 
non-zero of the divergence leads to that the state space 
volume is not conserved: Ascendancy is implied. The 
non-zero of the skew matrix, or the asymmetry of the se- 
lection matrix S, implies the existence of the transverse 
matrix T. 

Now, we give an explicit construction which demon- 
strates the existence and uniqueness connection between 
Eqs. (1,2) and Eqs. (8,9). Assuming that both Eq.(l) 
and (8) describe the same dynamics. The speed q t is 
then the same in both equations. The connection from 
Eq.(l) to (8) is straightforward: Multiplying both sides 
of Eq.(l) by [A(q t ,t) + T(q t , t)}- 1 leads to Eq.(8). Con- 
verting Eq.(8) into (1) is mathematically more involved. 

Using Eq.(8) to eliminate the speed in Eq.(l), and 
noticing that the dynamics of noise and the state variable 
behave independently, we have 

[A(q, t) + T(q, t)]f (q, t) = V^(q, t) , (11) 

and 

L4(q,i)+T(q,i)]C(q,i)=£(q,i) . (12) 

Here we have dropped the subscript t for the state vari- 
able, because time t is now a parameter. Those two equa- 
tions suggest a rotation in state space. 

Multiplying Eq.(12) by its own transpose of each side 
and carrying out the average over stochastic drive, we 
have 

[A(q,t) +T(q,t)]D(q,t)[A(q,t) -T(q,i)] = A(q,t) . 

(13) 

In obtaining Eq.(13) we have also used Eq.(2) and (9). 
Eq.(13) suggests a duality between the standard stochas- 
tic differential equations and Eq.(l): A large ascendant 
matrix implies a small diffusion matrix. It is a generaliza- 
tion of the Einstein relation in the case of zero transverse 
matrix 37 . 

Next we define an auxiliary matrix function 

G(q,t) = [A(q,t)+T(q,t)}- 1 . (14) 

Here the inversion '— 1' is with respect to the matrix. 
Using the property of the Wright fitness function ip: V x 
Vip = [(V x = (VjVj - ViVj)V> ], Eq.(ll) leads 

to 

V x [G-^q)] = , (15) 



which gives n(n — l)/2 conditions. The generalized Ein- 
stein relation, Eq.(13), leads to the following equation 

G + G T = 2D, (16) 

which readily determines the symmetric part of the aux- 
iliary matrix G, another n(n +1)/ conditions. The aux- 
iliary function may be formally solved as an iteration in 
gradient expansion: 

G = D + Q, (17) 

with Q = Hindoo AGj , AGj = 

EZii-^K^YDjS^ + iS^DjS 1 ], D = DS-S T D, 
Dj>x = (D + AGj-!) {[V x (D- 1 + AGj\)]f} (D - 
AGj_i). At each step of solving for AGj only linear 
algebraic equation is involved. One can verify that the 
matrix Q is anti-symmetric. For a simple case a formal 
solution of such algebraic equation was given in 38 , and 
an explicitly procedure was found for generic cases in 39 . 
Eq.(17) is a result of local approximation: If the selction 
matrix S, the diffusion matrix D are constant in space, 
the exact solution only contains the lowest order contri- 
bution in gradient expansion: Q = AGj = AGo- We 
regard Eq.(17) as the biological solution to Eq.(15) and 
(16), because it preserves all the fixed points of deter- 
ministic drive f. The connection from Eq.(8) to (1) is 
therefore uniquely determined: 

tf(q,t) = f c dq>-{G-\q>)i(q>)} 
A(q,t) = [G- 1 (q) + (G-)- 1 (q)]/2 . (18) 
T(q,t) = [G _1 (q) - (G r )" 1 (q)]/ 2 

Here the sufficient condition det(^4 + T) ^ is used, and 
the end and initial points of the integration contour G 
are q and qo respectively. 

We point out that in the absence of stochastic drive, 
i.e., e = in Eq.(2) and (9), above connection remains 
unchanged. 

B. Fokker-Planck Equation 

In many experimental studies in biology, a question is 
often asked on the distribution of the state variable as 
a function of time instead of focusing on the individual 
trajectory of the system. This implies that either there 
is an ensemble of identical systems or repetitive experi- 
ments are carried out. To describe this situation, we need 
a dynamical equation for the distribution function in the 
phase space. This goal can be accomplished by the so- 
called Fokker-Planck equation, or the difussion equation, 
or the Kolmogorov equation 14,15,22 . 

In this subsection, another procedure to find the equa- 
tion for distribution function is presented. It is natural 
from a theoretical physics point. This procedure will es- 
tablish that the Wright fitness function ip in Eq.(l) in- 
deed plays the role of potential energy in the manner 
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envisioned by Wright, and the steady state distribution 
will be indeed given by Eq.(3). Our starting point will 
be the second law, Eq.(l), not the standard stochastic 
differential equation, Eq.(8), from which most previous 
derivations started. 

The existence of both the deterministic and the 
stochastic drives in Eq.(l) suggests that there are two 
well separated time scales in the system: the microscopic 
or fine time scale to describe the stochastic drive and 
the macroscopic or course time scale to describe the sys- 
tem motion. The former time scale is much smaller than 
the latter. This separation of time scales further sug- 
gests that the macroscopic motion of the system has an 
" incrtial" : it cannot response instantaneously to the mi- 
croscopic motion. To capture this feature, we introduce a 
small constant inertial "mass" to and a kinetic momen- 
tum vector p for the system. Our state space is then 
enlarged: It is now a 2n-dimensional space. The dynam- 
ical equation for the system takes the form: 

q t = pt/m , (19) 

which defines the kinetic momentum, and 

Pt = -[A(q t ,t) + T{dt,t)]p t /m + W(qt, t) + £(qt, t) , 

(20) 

which is the extension of Eq.(l). We note that there is 
no dependent of ascendant matrix A and the stochastic 
drive on the kinetic momentum p. The Fokker-Planck 
equation in this enlarged state space can be immediately 
obtained 22 : 



I m 



• V + f • V r 



-+Vpl}p(q,p,t)=0. 

TO J J 

(21) 



Here f = pT/m + V q V>, and t, q, and p are independent 
variables. The subscripts in the d and V indicate the dif- 
ferentiation with respect to indicated variable only. The 
stationary distribution can be found, when the Wright 
fitness function is time-independent and bounded above, 



as 
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p(q, p, t 



oo) = — exp 



p 2 /2to - v(q) 



(22) 



with Z = J nlli <ki Uti dP* exp{-[p 2 /2TO - V(q)]/e} 
the partition function in the extended state space. There 
is an explicit separation of state variable and its kinetic 
momentum in Eq.(22). The elimination of the momen- 
tum in the small mass limit will not affect this dis- 
tribution. Hence, Eq.(22) confirms that the expected 
Boltzmann-Gibbs distribution, Eq.(3) from the Eq.(l) 
and (2), is the right choice. 

We proceed to outline the procedure to find the Fokker- 
Planck equation corresponding to Eq.(l) and (2) with- 
out the kinetic momentum p. We first illustrate how 
to recover Eq.(l) from Eq.(19) and (20). In the limit 



of m — > 0, the fast dynamics of kinetic momentum p t 
can always follows the motion of slow dynamics of state 
variable q t . Hence we may set p t = in Eq.(20) and re- 
place the kinetic momentum using Eq.(19), which is then 
Eq.(l) after moving the speed to the left-side of equation. 
For the Fokker-Planck equation, the explicit separation 
of the kinetic momentum and state variable in the sta- 
tionary distribution gives the guidance on the procedure: 
The resulting Fokker-Planck equation must be able to 
reproduce this feature. The Fokker-Planck equation is 
then found as 

d t p(d, t) = V T [-f (q) - Af (q) + D(q)V]p(q, t) , (23) 

with Af the solution of the equation V • Af + Af • V^i> — V • 
[GTG T Vip] = 0. If the probability current density is de- 
fined asj(q,t) = (f+Af-DV)p(q,i), the Fokker-Planck 
equation is a statement of the probability continuity: 



$p(q.*) + v-j(q,t) = o 



(24) 



The stationary state corresponds to the condition V • 
j(q, t — oo) = 0. One may verify that the stationary 
distribution p(q, t = oo) in Eq.(3) is indeed the time 
independent solution of the Fokker-Planck equation: The 
stationary probability current 

j(q, t = oo) = (GTG T + Af)VV(q) p(d, t = oo) , (25) 

and V • j(q, t = oo) = 0. 

The connection between the standard stochastic dif- 
ferential equation and Fokker-Planck equation has been 
under intensive study by biologists, physicists, chemists, 
mathematicians, and others over last 70 years 14 ' 15,17,22 . 
However, there exists an ambiguity for the generic non- 
linear situation 22,35 . We attribute this ambiguity to the 
asymptotic nature of the connection in which a procedure 
must be explicitly defined: Different procedures will in 
general lead to different results. Biologically, it is a state- 
ment on how the dichotomy of deterministic and stochas- 
tic drives is done, a genuine indication of the hierarchical 
nature of the dynamics. What has been demonstrated in 
this subsection is one way of carrying out this procedure. 



C. Detailed Balance Condition 

There is an important class of evolution dynamics in 
which the anti-symmetric matrix Q = 0. Under this 
condition, the transverse matrix T = 0, and Af = 0. 
The Fokker-Planck equation becomes 



d t p(<l,t) = V T [-f(q)) + Z>(q)V]p(q,t) 



(26) 



and the stationary probability current is everywhere zero 
in state phase: 



j(q,i = oo) = . 
In this situation one may find that 



(27) 
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VV'(q) = J D- 1 (q)f(q) 



(28) 



and A = D^ 1 . The Wright fitness function and the 
connection between Eq.(l) and the standard stochastic 
differential equation can be directly read out from equa- 
tions. This is the well-known symmetric dynamics in 
biological 14,15,17 and physical 22 sciences. A condition to 
generate this kind of equilibrium state in biology was first 
noticed by Hardy 40 and Weinberg 41 . This zero probabil- 
ity current condition is usually called the detailed balance 
condition. 



IV. EXAMPLES 

In this section we discuss four examples. The first one 
is of predator-prey model like. In this model we illus- 
trate how to approximately compute the Wright fitness 
function ip, the ascendant matrix A and the transverse 
matrix T, that is, how to make the connection between 
the conventional formulation and Eq.(l) and (2). Second 
and third examples are current models in evolutionary 
biology. Fourth example is the usual diffusion equation. 
We use them to demonstrate that they may be discussed 
within the present mathematical formulation. 



A. Predator-Prey Model 

The example whose dynamical equation is in the form 
of standard stochastic differential equation is the generic 
predator-prey process. Under the diffusion approxima- 
tion, both the diffusion matrix D and the deterministic 
drive can be obtained from the master equation. The 
diffusion approximation is valid when a large number of 
birth and death events occurs on the macroscopic time 
scale 22 . 

We remark that in reality, both the intrinsic and the 
extrinsic noise coexist. They are equally important and 
measurable 22 . The stochasticity- ascendancy relation, or 
the fundamental theorem of evolution, treats both of 
them on the equal footing to determine the connection 
between Eq.(l) and (8). 

We now give an explicit demonstration of how to ob- 
tain Eq.(l) from Eq.(8). Here q\ and q 2 represent num- 
bers of two species in a habitat. We assume the spatial 
distribution is uniform. The deterministic drive f con- 
sists of two positive terms, birth and death: 



/»(q) = /»b(q) - /»d(q) * = i, 2 



(29) 



with the subscripts b and d stand for the birth and 
death respectively. Under the diffusion approximation, 
the stochastic drive i s 14 > 42 > 22 

Ci(q,t) = ^ fi P (ci)(i P (t) + V/w(q)CwW * = 1,2, (30) 

with Cip(t)i Qd(t) are unity random variables and possible 
correlation among them. Therefore the diffusion matrix 



D can be readily obtained, which is what needed below. 
We remark that the equation similar to predator-prey 
equation has been emerged in the study of the robustness 
of the gene regulatory network of phage A 43 . 

The construction of Eq.(l) from Eq.(8) will be given to 
the lowest order in the gradient expansion. The useful- 
ness of this approximated construction can be illustrated 
for following two reasons. First, in many practical appli- 
cations, lowest order approximation is already enough 43 , 
because it is exact in the strictly linear case. Second, sev- 
eral salient features of the connection becomes apparent 
without undue mathematical complications. An impor- 
tant quantity is the selection matrix S. According to the 
definition following Eq.(13), 

5*11 = V1/1 , S12 = V2/1 , S21 = V1/2 , S 22 — V2/2 • 

(31) 

Eq.(10) will not change under the gradient approxima- 
tion. In the lowest order gradient approximation, Eq.(9) 
becomes simple. We collect them here: 



GS T -SG T = 
G + G T = D 



(32) 



In two dimensions the matrix manipulation is particu- 
larly straightforward. We note that any 2x2 matrix M 
can uniquely decomposed in terms of Pauli matrices, <7j 
with i = 1, 2, 3, and the identical matrix 1: 

M = M\(T\ + M 2 a 2 + M30-3 + tr(M)/2 1 , 



with tr denotes the trace and o\ 
1 



, and here i 




-i 

1 y ai ~ yo -i 
this relationship, the equation for antisymmetric part of 
the auxiliary matrix G = D + Q from Eq.(32) is 



QS T + SQ = (SD - DS T ) 



(33) 



Using the matrix decomposition and the properties of 
Pauli matrices, we obtain 



Q= (SD- DS T )/tr(S) . 
Note that for the 2x2 matrix M 



(34) 



Mn M12 
M 2 i M 22 



_L / M 22 -M12 

det(M) ^ -M21 Mn 



The ascendant matrix A and the transverse matrix T can 
be found according to Eq.(8): 



' tf(q) 
Ml) 

k T(q) 



Ldq'-IG^qW)] 
[ D 22 -D 12 
\-D 12 D n 
-Q/det(G) 



/det(G) 



(35) 
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In two dimensions, det(G) = det(D) + det(Q) = 
D22D11 — D\ 2 + Qi 2 and is obviously non-negative. 

To summarize, the computation of quantities in Eq.(l) 
from Eq.(8) is as follows: First, to establish Eq.(8) from 
the biological problem, continuous approximation should 
be used. If the problem is given in terms of master equa- 
tion, the usual diffusion approximation will be employed, 
which gives both the diffusion matrix and the determinis- 
tic drive 15,22 . After this is done, the procedure prescribed 
in section III. A is employed to find the ascendant matrix, 
the transverse matrix and the fitness matrix. According 
to the biological problem, an additional approximation 
may be used to reduce computation effort but with the 
desired accuracy, as done here as well as in Zhu et al. 43 . 

B. Bateson-Dobzhansky-Muller Model 

Based on the verbal descriptions of Bateson, Dobzhan- 
sky, and Muller, Gavrilets has developed a mathemat- 
ical description of adaptive landscape with many loci 
and traits to describe the speciation: the holey adap- 
tive landscapes 11 . It puts the emphasis on narrow fitness 
bands linking large fitness areas with the same fitness. 
The fitness peaks, expressed at rugged adaptive land- 
scapes, are argued to be not important when number of 
traits and loci is large, contrast to the demonstration by 
Kisdi and Geritz 12 . 

This is a very attractive idea, and is consistent with 
present framework. With an appropriate coarse grain av- 
erage, it may be possible to compute this holey adaptive 
landscape in terms of the Wright fitness function defined 
in the present article. The speciation may be interpreted 
as the diffusion along the narrow fitness bands, a process 
may be much faster than that described by Eq.(5) for the 
peak hopping. 

C. Symmetry-Breaking Model 

The symmetric evolution dynamics was explicitly dis- 
cussed by Stewart for speciation 10 : 

q t - V0(q t ) . (36) 

This is a deterministic equation. Stochastic model was 
also mentioned by Stewart 10 . 

In the light of present discussion, the diffusion ma- 
trix D is equivalent to the identity matrix 1. This is 
a case satisfying the detailed balance condition. Hence 
the Wright fitness function can be directly obtained: 
V'(qt) = <Hqt)> an d the ascendant matrix A = 1, and 
transverse matrix T = 0. The steady state distribution 
is then given by the Boltzmann-Gibbs like distribution, 
Eq.(3). All the statistical physics methodology can then 
be applied here. Naturally one may make use of the 
idea of symmetry-breaking as a way for self-organization 



needed for speciation. We refer readers to the beautiful 
discussion presented by Stewart 10 . 

We remark that one may even make use of the self- 
consistent mean-field approximation, a powerful math- 
ematical tool in statistical physics 44 , to search for the 
indication of symmetry-breaking. 

D. Diffusion Equations 

Expressing evolutionary dynamics in the form of 
diffusion equations has been common in population 
genetics 15,17 . Particularly, the one-dimensional diffusion 
equation has been thoroughly studies. It belongs to the 
class satisfying the detailed balance condition. Hence, 
as discussed in subsection III.C, the transverse matrix 
is zero, and the ascendant matrix and the gradient of 
the Wright fitness function, or adaptive function, can be 
readily identified. The results, such as the equilibrium 
population, are consistent with what known. We remark 
that if there is an overall factor difference in terms of 
the diffusion matrix, one needs to examine which di- 
chotomy of deterministic and stochastic drives is used: 
Ito, Stratonovich, the present prescription, or something 
else. A discussion of this feature can be found in the 
monograph of van Kampcn 22 . 

In case the detailed balance condition would not hold 
for higher dimensional diffusion equations, a way to pro- 
ceed along the present formulation is discussed at the end 
of subsection IV. A. 



V. DISCUSSIONS 

Before further going to further discussion on the impli- 
cation of the present mathematical formulation, it should 
be kept in mind that above three laws must be regarded 
as mathematically what the evolutionary dynamics might 
be. They are by not means the exact description. Hav- 
ing made this statement, we nevertheless remark that 
although proposed three laws for evolutionary dynamics 
are based on the continuous approximation in terms of 
time and population, it is possible that main features 
discussed in the present article may survive in discrete 
cases. 

In the present attempt to unify approaches from biol- 
ogy and physical sciences, the possible existence of laws 
such as expressed by Eq.(l) and (2) should not be too 
surprised if one views them from two important prin- 
ciples which have been rigorously validated: 1) Simple 
equations can generate extremely complicated patterns 
and phenomena 31 ' 30 ' 29 ; 2) Each level of description has 
its own laws which cannot be derived in a naive reduc- 
tionism manner 45 ' 44 ' 46 . The connections between levels 
are asymptotic and emerging phenomena frequently oc- 
cur at higher levels. 
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There are several quantitative advantages in present 
formulation of evolutionary theory. With the Wright fit- 
ness function defined as in Eq.(l), an independent way to 
calculate it is obtained, as one may follow what indicated 
in Section III. This adds more predictive power into the 
evolutionary dynamics. 

As expressed by Eq.(2) and discussed in section 
2.4, Fisher's fundamental theorem of natural selection 
(Fisher's theorem) becomes transparent and indispens- 
able. This may provide a much-needed step to better 
understand Fisher's great insight. 

Combining both Fisher's and Wright's insights, the 
Wright fitness function and Fisher's theorem provide a 
quantitative measure to discuss robustness, stability, and 
the speciation. Eq.(5) is such an example. 

It is interesting to point out the remarkable similar- 
ity between the adaptive landscape of Wright 19 and the 
developmental landscape of Waddington 47 . The present 
mathematical formulation can deal with both cases. An 
example on the gene regulatory network in phage A has 
already been studied by Zhu et al. 43 . This suggests a 
unification between genetics and developmental biology. 

However, one may ask why to use Eq.(l) and (2) in- 
stead of more conventional Eq.(8) and (9): After all their 
equivalence has been demonstrated above. Here we offer 
three reasons to favor Eq.(l) and (2): 

i) Quantities presented in Eq.(l) can be directly related 
to experimental observation. For example, Eq.(3) gives 
a direct connection between the Wright fitness function 
and the population density in steady state. By observing 
the dynamical behaviors, information on the ascendant 
and transverse matrices can be obtained. Also, Eq.(5) 
can relate stability to the Wright fitness function. This 
direct contact with experimental data is an indication of 
the autonomy of the focused level of description. 

ii) Eq.(8) and (9) lack the visualizing ability for the 
global dynamics behavior. For example, in a nonlinear 
dynamics with multiple local maxima, it is not clear from 
Eq.(8) and (9) which maximum is the largest one, and 
how easy it might be to mover from one maximum to 
another. One could find this answer by a direct real time 
calculation. But this is usually computationally demand- 
ing, if not impossible. 

iii) Eq.(l) and (2) give an alternative modeling of evo- 
lutionary dynamics, which can be advantageous in cer- 
tain situations. For example, the direct use of fitness 
function in Stewart's modeling 10 makes the symmetry- 
breaking idea very transparent from statistical physics' 
point of view. 

Finally, we point out an interesting mutual reduction 
loop between biology and physics. In the discussion of 
the first law we have remarked that one may regard the 
Newtonian dynamics as a special case of the first law. 
This implies that the Newtonian conservative dynamics 
is a special case of the present second law, hence the 
Darwinian dynamics. Here the opposite statement also 
exists: Under an appropriate condition equations in the 
form of Eq.(l) and (2) can be derived from the Newto- 



nian dynamics 50 , therefore the Darwinian dynamics 
may also be regarded as a special case of the Newtonian 
dynamics. 

VI. CONCLUSIONS 

Based on continuous approximation we have postu- 
lated three laws to mathematically describe the evo- 
lutionary dynamics. The most fundamental equation, 
the second law, has been expressed in a unique form 
of stochastic differential equation. Four dynamical ele- 
ments have used in our formulation: the ascendant ma- 
trix, the transverse matrix, the Wright fitness function, 
and the stochastic drive. We have demonstrated that 
present laws are consistent with previous approaches in 
biology, but appear more suitable to discuss stability and 
other phenomena quantitatively. Various important re- 
sults, such as Fisher's fundamental theorem of natural 
selection and Wright's adaptive landscape, as well as the 
developmental landscape, are apparently unified in the 
present formulation. It appears that the present quanti- 
tative formulation has captured the main aspects of the 
Darwinian dynamics. 
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